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ABSTRACT 

We present three-dimensional numerical simulations of particle clumping and planetesimal formation in 
protoplanetary disks with varying amounts of solid material. As centimeter- size pebbles settle to the mid- 
plane, turbulence develops through vertical shearing and streaming instabilities. We find that when the pebble- 
to-gas column density ratio is 0.01, corresponding roughly to solar metallicity, clumping is weak, so the pebble 
density rarely exceeds the gas density. Doubling the column density ratio leads to a dramatic increase in 
clumping, with characteristic particle densities more than ten times the gas density and maximum densities 
reaching several thousand times the gas density. This is consistent with unstratified simulations of the streaming 
instability that show strong clumping in particle dominated flows. The clumps readily contract gravitationally 
into interacting planetesimals of order 100 km in radius. Our results suggest that the correlation between host 
star metallicity and exoplanets may reflect the early stages of planet formation. We further speculate that 
initially low metallicity disks can be particle enriched during the gas dispersal phase, leading to a late burst of 
planetesimal formation. 

Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary disks — 
solar system: formation — turbulence 



l. introduction 

The concentration of particles to high spatial densities pro- 
motes the formation of planetesimals, the super-kilometer 
scale building blocks of planets. Drag forces on pebbles 
and rocks in disks lead to spontaneous particle clumping 
dGoodman & Pindortl2000l). The discovery of a linear stream- 
ing instability (I Youdin & Goodm an 2005) shows that clump- 
ing is a robust consequence of particles drifting in and 
gas flowing out i n disk s with some radial pre s sure s upport 
dNakagawa eTafl 119861) . Uohansen & Youdinl (120071) stud- 
ied the non-linear saturation of the streaming instability, ne- 
glecting vertical gravity and self-gravity. Those simulations 
showed that groups of boulders accelerate the gas around 
them towards the Keplerian velocity, reducing the radial drift 
locally a nd leading to tempor ary concentrations of boulders 
(s ee also 1 Johansen et al. 20 061). 

IChiangl (|2008|) and iBarrancol (|2009|) recently performed 
three-dimensional (3D) simulations of vertical shear instabili- 
ties in Keplerian disks in the single fluid limit where particles 
and gas have exactly the same velocities. These studies con- 
firmed expectations that mid-plane turbulence develops when 
the Richardson number Ri < 1 . While perfect coupling is a 
good approximation for small grains, it cannot include verti- 
cal settling or in-plane streaming motions. 

In this Letter we present 3D simulations of the motion of 
gas and pebbles in sub-Keplerian disks, including vertical 
gravity and particle sedimentation. Thus, we can study the 
combined effect of vertical shearing and streaming instabil- 
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ities, as particles self-consistently settle towards - and are 
stirred from - the mid-plane. We exclude external sources of 
turbulence, including magnetorotation al instabilities ( which 
can actually prom ote clumping, see Joha nsen et al.l 120071: 
iBalsara et al.l 120091) . Our hydrodynamical simulations of- 
fer a first approxima tion to dead zones with low ioniza- 
tion (ISano et al.l |2000|) where turbu lent surface layers drive 
only weak motio ns in the mid-plane ( Fle ming & Stone 1120031 : 
lOishi et al.ll2007l) . 

In this non-magnetized limit, we in vestigate the clumpin g 
of smaller particles than considered in I Johansen et al.l (120071) . 
which increases the likelihood of coagulation up to the ini- 
tial sizes. We find that clumping of pebbles in the mid-plane 
layer increases sharply above a threshold mass fraction of 
solids roughly consistent with solar metallicity. Thus plan- 
etesimal formation may help explain the high probability of 
finding gia nt exo planets aroun d star s rich in heavy ele ments 
(lGonzalezlll997l:jSantos et al.l l2001l: iFischer & Valentil 120051: 
Uohnson & Appl2009l) . 

2. SIMULATIONS 

We perform 3D hybrid simulations. They model gas on a 
fixed grid and solids with superparticles, each representing 
a swarm of actual particles. We solve the standard shearing 
sheet dynamical equations for a frame rotating at the Keple- 
rian frequency £1 at a fixed orbital distance r from the star. 
The axes are oriented such that x points radially outwards, y 
points in the orbital direction, while z points vertically out of 
the disk. The gas is subject to a global radial pressure gradient 
that reduces the gas orbital speed by Av = -0.05c s w 25 m s -1 
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Fig. 1. — Results of 128 3 sedimentation simulations for various values of the pebble abundance Z p (different columns). The upper row of images shows the 
particle density averaged over the azimuthal direction at the end of each simulation. The wave-like pattern in the particle density is a consequence of coherent 
radial drift and vertical oscillations. The lower row shows particle density averaged over both y and z directions, as a function of radial coordinate x and time t. 
At roughly solar metallicity (left column) the particle surface density is uniform. Super-solar metallicities (middle and right columns) produce significant clumps 
that merge over time because weaker clumps drift faster. 



(IChiang & Goldreichlll997l) . The sound speed c s , gas scale 
height H g = c s /£l and mid-plane gas density po are the natural 
units of the simulation. 

The motion of gas and particles are coupled through 
momentum-conserving drag forces with particle friction time 
Tf. Our dynamic a l equ ations are identical to those of 
lYoudin & Johansen! (120071) . with the addition of a vertical 
gravitational acceleration g z = -Q?z affecting both gas and 
particles. 

The superparticles are evenly distributed in mass and 
number into four bins of normalized friction time £lrf = 
0.1,0.2,0.3,0.4. These friction times are characteristic of 
compact solids with radius a p w 3, 6,9, 12 cm at r = 5 AU in 
the Minimum Mass Solar Nebula (Weidenschillin j [T977allbl: 
lHavashil[T981l) . Rescaling to r = 10 AU yields a v w 1-4 cm. 
We colloquially refer to this range of particle sizes as pebbles 
to contrast with larger Clrf ^1.0 boulders. The total pebble 
mass is fixed by setting the pebble-to-gas column density ra- 
tio Z p = (Zp)/ (Z g ), where (Z p ) and (Z g ) are the mean particle 
and gas column densities, taking into account that most of the 
gas resides beyond the vertical extent of the simulation box. 
This pebble abundance turns out to be the crucial parameter 
for triggering particle clumping. 

The total abundance of condensable m aterials beyond the 
ice line was estimated by iHayashl (1198 lb to be Z c w 0.018, 
while more up-to-date models give a somew hat lower value 
of Z c g 0.015 at temperatu res less than 41 K (lLoddersll2003l: 
lAUende Prieto et al.ll200lh . For our models a greater uncer- 
tainty is the efficiency of conversion from dust grains to peb- 
bles. Assuming that a majority 2/3) of the condensable 
solids are bound in pebbles, Z p = 0.01 corresponds to so- 
lar metallicity. We also experiment with higher values of 
Z p = 0.02 and Z p = 0.03, which are motived both by stars with 



super- solar metallicities and by mechanisms that enrich the 
solids-to-gas ratio in disks (see A given pebble abun- 
dance would correspond to higher values of the metallicity if 
pebbles make up a smaller fraction of the condensable mate- 
rial. 

We use a box size of L x = L y = L z = 0.2H g and resolutions 
of 64 3 zones with 125,000 particles, and 128 3 zones with 
1,000,000 particles. This relatively small box size is cho- 
sen to capture typical wavelengths of streaming and Kelvin- 
Helmholtz instabilities. The gas density is in vertical hydro- 
static equilibrium. Particle positions are initialized to give 
a Gaussian density distribution around the mid-plane with 
scale height H p = 0.02// g , while gas and particle velocities 
ar e initially set to match the drag force equilibrium solution 
of lNakagawa et alJ (119861) . 

3. PARTICLE CLUMPING 

Since the disk is initially laminar, particles settle to the disk 
mid-plane. As particles collect in the mid-plane, they accel- 
erate gas there towards the Keplerian orbital speed. This gen- 
erates vertical shear that can drive Kelvin-Helmholtz insta- 
bilities. The velocity difference between gas and solids also 
triggers streaming instabilities. The resulting turbulence halts 
particle sedimentation and can lead to clumping. 

In Fig. [T] we show results for pebble abundances Z p = 
0.01 (roughly solar), Z p = 0.02 (super-solar) and Z p = 0.03 
(strongly super-solar). The Z p = 0.01 simulation shows min- 
imal clumping, but reveals instead a surprising interaction 
between settling and stirring. The particle layer is not cen- 
tered uniformly on the mid-plane, as usually assumed. In- 
stead the solids organize in a wave-like pattern (top left panel 
of Fig. [T]) produced by coherent vertical oscillations com- 
bined with radial drift. We note that the streaming instabil- 
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Fig. 2. — Cumulative distribution functions of the relative speed between particles of various sizes, weighted by the relative speed itself. Vertical lines mark 
50% and 10% in the distribution functions. The typical collision speed lies between 0.5 and 5 ms -1 , although collisions at up to 10m s _1 occur as well. Collision 
combinations 0.2+0.2, 0.2+0.3, 0.2+0.4 showed very similar results and are not included in the figure. The low resolution simulation (dotted lines) shows a bit 
higher collision speeds, because less clumping occurs in it. The lower left panel shows the mean collision speed as a function of the density of particles in the 
grid cell. Grid cells with a high particle density have typical collision speeds below 2 ms" 1 , At low particle densities the collision speed between Q.Tf = 0.1 and 
Qrf = 0.4 particles rises, in agreement with differential drift (NSH analytical solution, Nakagawa et al. 1986), while the overall collision speed tends towards the 
collision speed between the Q.Tf = 0.1 bodies that dominate low density regions. 



ity is kno wn to drive coherent verti cal motions in clumps of 
boulders dJohansen & Youdinl 120071) . The characteristic par- 
ticle density p p = {p p ) / {p p ) ~ 0.9p g is high enough to exert 
feedback on the gas. However the average mid-plane density 
pjmid) - o.6p g is diluted by the voids in the wave-like pattern, 
and is too small to trigger strong radial clumping. 

Increasing the pebble abundance to Zp = 0.02 and Z p = 0.03 
triggers strong overdensities in the particle layer (see Fig.Q]). 
Initially, streaming instabilities produce many azimuthally 
extended clumps in the simulation box. Since the denser 
clumps drift inward more slowly, mergers result in a sin- 
gle dominant clump. The characteristic particle density in- 
creases due to the strong clumping, with p p /p g ~ 0.9, 18,73 
for Z p = 0.01,0.02,0.03, respectively. The maximum density 
in the box increases from less than ten to more than 2,000 
times the gas density when the pebble abundance is increased 
from Z p = 0.01 to 0.02. 

So why does the clumping increase so sharply for order 
unity changes to the pebble abundance Z p ? Fig. 0] (top pan- 
els) shows that the vertical extent of the wave- shaped par- 
ticle layer decreases with Z p , because less kinetic energy 
is rel eased in particle-dominated flows dJohansen & Youdinl 
120071) . Consequently the average mid-plane particle 
density increases, with Pp^/Pg = 0.6,2.0,9.0 for Z p = 
0.01,0.02,0.03. Strong clumping is thus expected for the 
higher metallicity cases, because the streaming instability is 
more powerful — for linear growt h and nonlinear clumping— 
when the average P v /pe > 1 (lYoudin & GoodmanT 120051: 
iJohansen & Yo udin 2007). A runaway develops as clump- 
ing limits stirring, giving higher average densities and in turn 
more clumping. 

We can estimate the transitional value Z* above which par- 
ticles strongly clump. Assume that the sub-layer thickness is 



set by vertical shear turbulence to beH p ~ Av/(2Q) as appro - 
priate for small particles (ISekivall 199 a lYoudin & Shull2002l) . 
Then the clumping threshold pjf 11 ^ > p g requires 



For Av = 0.05c s this yields a critical pebble abundance of 
Z* ~ 0.025. The slight overestimate of Z* may arise be- 
cause our pebbles settle more effectively than small grains. 
The result that clumping happens for lower Z p when Av/c s is 
smaller ( weaker radial pressure support) is seen in the simu- 
lations of IJohansen et al.1 (120071) . 

The same metallicity thres hold arises in p article layers with 
constant Richardson number (ISekiyalll998l) . Above Z* a high 
density cusp forms in the mid-plane because the shear tur- 
bulen ce cannot lift solids w ith mass exceeding the local gas 
mass (lYoudin & Sh u 2002). Vertical self-gravity causes the 
analytic cusps to diverge, while streaming instabilities are re- 
sponsible for the strong clumping in our dynamical simula- 
tions. 

4. COLLISION SPEEDS 

The growth of solids into planetesimals depends on colli- 
sion spe eds, which affect the r ates of fragmentation and coag- 
ulation (iBlum & W urm 2008). The threshold speed for catas- 
trophic destruction of a meter- scale body, with the materia l 
strength of compact rock, is around 10 m s -1 (lBenzll2000l) . 
Porous aggregates, with less material strength, have an even 
lower threshold for destruction. Sticking of equal- sized bod- 
ies > 1 mm is not observed. However, small grains can adhere 
at speeds less than ~lms _1 (Blum & Wurm 2000). 

From our simulation snapshots we measured relative speeds 
of all particle pairs within the same grid cell, giving 
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(1/2)N(N- 1) unique collisions in a grid cell with N parti- 
cles. All collisions are assumed to be head on; for a random 
distribution of impact parameters the normal velocity would 
be reduced by 50% compared to the relative speed. We av- 
eraged over snapshots taken between 30 and 50 orbits of the 
64 3 and 128 3 runs with Z p = 0.02. 

The distribution of collision speeds between particles of 
various sizes is shown in Fig. [2] We weight the occurrence 
of relative speeds by the speed itself to reflect that high speed 
collisions take place more frequently than low speed colli- 
sions. The typical collision speed lies in the interval 0.5- 
5 ms -1 for a sound speed c s « 500 ms -1 . Fig. [2] also shows 
that the collision speed is systematically lower in the denser 
particle clumps where the typical speed falls to less than 
2ms -1 . Averaging collision speeds without weighting by the 
speed itself leads to much lower typical speeds, around 0.2- 
2 ms -1 . 

The low collision speeds are surprising given the sub- 
Keplerian speed of 25 m s -1 . The differential radial drift alone 
is 12.7 ms -1 between Clrf = 0.1 and Clrf = 0.4 test particles. 
However, the large bodies reside in dense clumps in the mid- 
plane layer where the differential radial drift is much reduced 
and the turbulence is weakened by mass loading (see insert in 
Fig. E]). 

In simulations of particles in magnetorotational turbulence, 
collision speeds are measured to be much higher, of order 
^lOms -1 (and well beyond for different- siz ed particles, see 
I Johansen et all 120071: iCarballido etaH 12008b . The more be- 
nign collision environment in the current simulations is more 
resistant to destruction and may even allow for growth by 
sticking. Low random velocities in high density clumps also 
favor gravitational collapse on a much shorter time-scale. 

5. SELF-GRAVITY 

Collective gravi ty can act to assemble solid particles 
into planetesimals (Safronov 1969; Goldreich & Ward 19731; 
lYoudin & ShJl2002HJohansen et alJl2007t iCuzzi et alJl2008h . 
The clumps in our Z p = 0.02 and Z p = 0.03 simulations reach 
densities beyond the Roche density, 

<* = |^F~ 100 *' (2) 

so self-gravity will be dynamically relevant. We activate the 
self-gravity in a Zp = 0.02 simulation that already developed 
for ~40 orbits, starting at a chosen time that displayed rel- 
atively weak clumping. A few orbits later a dense filament 
forms due to the streaming instability, but the concurrent ac- 
tion of the self-gravity leads to the condensation of seven 
gravitationally bound objects (see Fig. 0. Initially the plan- 
etesimals actively accrete pebbles from their surroundings. 
However, the accretion slows down after 5 orbits as 20% of 
the solid mass is locked in seven clumps with masses corre- 
sponding to solid bodies of radii between 100 and 200 km. 
Further growth can occur both by accreting dust and pebbles 
and by planetesimal-planetesimal impacts, since the modest 
gas density fluctuations in mid-plane layer turbulence can not 
pump c ollision speeds of large planetesi mals into the erosive 
regime (llda et al.ll2008l: I Yang et al.ll2009l) . 

6. CONCLUSIONS 

In this Letter we present numerical evidence that the 
streaming instability, which arises in the coupled motion of 
gas and solids, can alone provide the necessary ingredients 
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Fig. 3. — Snapshots of particle column density after self-gravity is acti- 
vated. A dense nearly axisymmetric filament forms by the concurrent ac- 
tion of streaming instability and self-gravity, but it continues to fragment into 
gravitationally bound clusters of pebbles and rocks. The clusters initially 
show clear signs of active accretion of solids from the surroundings (top in- 
sert). After 5 orbits approximately 20% of the solid mass is in seven bound 
clusters, each containing a mass equivalent to a 100-200-km scale solid body. 

for the formation of planetesimals. With a high enough abun- 
dance of pebble- sized objects, the streaming instability drives 
the formation of high-density clumps in the sedimented mid- 
plane layer. Collision speeds inside those clumps are rela- 
tively low, around a few meters per second. The clumps get 
dense enough to contract gravitationally on a short time- scale 
into gravitationally bound clusters containing mass equivalent 
to solid planetesimals with radii of order 100 km. 

Our simulations display a threshold dependence of clump- 
ing on the global value of the pebble abundance. At a peb- 
ble abundance of Z p = 0.01, corresponding roughly to solar 
metallicity, powerful turbulence develops and puffs up the 
mid-plane layer to such widths that clumping is strongly sup- 
pressed. This environment is not conducive to planetesimal 
growth, with low particle densities and high collision speeds. 
However, a doubling of the pebble abundance to Z p = 0.02 
leads to strong clumping with local peak particle densities ex- 
ceeding a thousand times the gas density. This threshold be- 
havior arises because the streaming instability triggers strong 
clumping when the particle density exceeds the gas density in 
the mid-plane. 

If the fraction of condensables contained in pebbles is 
less than we assumed in this discussion (approximately 2/3), 
then the threshold we describe corresponds to higher val- 
ues of the bulk metallicity. Coagulation-fragmentation mod- 
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els dBrauer et a l. 2008) and observ ations of pebbles in proto- 
planetary discs (IWilner et al.ll20Q5h will be crucial for putting 
tighter constraints on the pebble abundance and the location 
of the metallicity switch. 

The pebble abundance may increase as solid particles 
drift radially inward, convergi ng in the inner few dozen 
AUs of the protop lanetary disk (Stepinski & Vala geaslll996t 
lYoudin & Chia ng 2004). A higher Zp also will occur if the 
gas in the disk preferentially evaporates. For example, pho- 
toevaporation removes gas from the disk surface, leaving 
behind larger solids that have sedimented to the mi d-plane 
dThroop & Ballvl 120051: IXlexander & Armitagd 120071) . This 
opens a window for late planet formation with planetesimals 
forming as the gas dissipates. The planets resulting from plan- 
etesimals formed so late will not have time enough to accrete 
a large gaseous envelope before the ga s dissipates and woul d 
also suffer less severe type I migration (lYoudin & Shull2002|) . 
This could explain why super-Earths and exo-Neptunes seem 
able to form equally well around stars of low and high metal- 
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The correlation between host star metallicity and giant ex- 
oplanets can be partly attributed to the efficienc y of form- 
ing seyeral-Earth-mass cores from planetesimals (llda & Linl 
120041: iMordasini et alJ [20091) . Our findings show that high 
metallicity already facilitates planetesimal formation, giving 
further su pport to the core ac cretion model for giant planet 
formation dPollack et al.ll 19961) . 
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